This is the 1D diffusive testcase against observed soil moisture dynamics and two Richards solvers of the particle model random walk engine

In [2]:

    

import numpy as np
import pandas as pd


    

In [3]:

    

obs=pd.read_csv('soil_moisture_weiherbach.csv')
obs.index=pd.to_datetime(obs['DATE'] + ' ' + obs['TIME'],format='%d/%m/%y %H:%M')
obs = obs.drop(['DATE','TIME'], 1)


    

In [4]:

    

#define plot function for soil moisture data
def hydroplot(obs,mlab,mlabcols,fsize=(6, 6),saving=False,upbound=40,lowbound=10,catch='Catchment',precscale=100.,cumprec=False,align=False,tspacehr=6,ccontrast=False,descriptor='Sensor\nLocation'):
    '''
    This is a rather simple function to plot hydrological data (soil moisture and precipitation) of a pandas data frame.
    It is based on some excellent examples by Randy Olson and may need heavy adaption to your data.
    (CC BY-NC-SA 4.0) jackisch@kit.edu
    
    fsize: Provide figure size as tuple.
    saving: Provide a file name if you want it saving.
    XXbound: Give bounds of left axis.
    catch: Provide catchment name.
    precscale: Scaling if your prec data is not in mm.
    cumprec: True: cumulative precipitation is plotted.
    
    The functions assumes a pandas data frame with a time stamp as row names.
    You may prepare this as:
    obs=pd.read_csv('soil_moisture_file.csv')
    obs.index=pd.to_datetime(obs['DATE'] + ' ' + obs['TIME'],format='%d/%m/%y %H:%M')
    obs = obs.drop(['DATE','TIME'], 1)
    Moreover, precipitation should reside in column 'Prec'
    '''
    
    # These are the "Tableau 20" colors as RGB.  
    tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),  
             (44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),  
             (148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),  
             (227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),  
             (188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]  
    # Scale the RGB values to the [0, 1] range, which is the format matplotlib accepts.  
    for i in range(len(tableau20)):  
        r, g, b = tableau20[i]  
        tableau20[i] = (r / 255., g / 255., b / 255.)  
        
    
    figure(figsize=fsize)  
      
    # Remove the plot frame lines. They are unnecessary chartjunk.  
    ax = subplot(111)  
    ax.spines["top"].set_visible(False)  
    ax.spines["bottom"].set_visible(False)  
    ax.spines["right"].set_visible(False)  
    ax.spines["left"].set_visible(False)  
      
    # Ensure that the axis ticks only show up on the bottom and left of the plot.  
    # Ticks on the right and top of the plot are generally unnecessary chartjunk.  
    ax.get_xaxis().tick_bottom()  
    ax.get_yaxis().tick_left()  
      
    # Limit the range of the plot to only where the data is.  
    # Avoid unnecessary whitespace.  
    ylim(lowbound, upbound)  
    
          
    # Make sure your axis ticks are large enough to be easily read.  
    # You don't want your viewers squinting to read your plot.  
    yticks(range(lowbound, upbound, 10), [str(x) + "%" for x in range(lowbound, upbound, 10)], fontsize=14)  
    xticks(fontsize=14)  
    
    ax.xaxis.set_minor_locator(matplotlib.dates.HourLocator(interval=tspacehr))
    ax.xaxis.set_minor_formatter(matplotlib.dates.DateFormatter('\n%H'))
    ax.xaxis.grid(False, which="minor")
    ax.xaxis.grid(True, which="major")
    ax.xaxis.set_major_locator(matplotlib.dates.DayLocator())
    ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('\n\n%d.%m.%y'))
    
    # Now that the plot is prepared, it's time to actually plot the data!  
    # Note that I plotted the majors in order of the highest % in the final year.  
    majors = mlab  
    colnames=mlabcols
    
    #get positions
    yposlab=obs[colnames].values[-1,:]
    yposlab=np.round(yposlab)
    if align:
        while any(np.diff(yposlab)==0.):
            overlap=np.diff(yposlab)
            overlapavoid=0.*overlap
            overlapavoid[overlap==0.]=1.
            overlap[np.where(overlap==0.)[0]+1][overlap[np.where(overlap==0.)[0]+1]==1.]=overlapavoid[np.where(    overlap==0.)[0]+1]+1.
            yposlab[1:]=yposlab[1:]+overlapavoid
    else:
        if any(np.diff(yposlab)<2.):
            yposlab=upbound-(np.arange(len(yposlab))*2.+15.)
    
    for rank, column in enumerate(majors):  
        # Plot each line separately with its own color, using the Tableau 20  
        # color set in order.
        if ccontrast:
            plot(obs.index, obs[colnames[rank]].values, lw=2.5, color=tableau20[rank*2+2])  
        else:
            plot(obs.index, obs[colnames[rank]].values, lw=2.5, color=tableau20[rank+2])  
    
        # Add a text label to the right end of every line. Most of the code below  
        # is adding specific offsets y position because some labels overlapped.  
        if ccontrast:
            text(obs.index[-1]+ datetime.timedelta(hours=0.1), yposlab[rank], column, fontsize=14, color=tableau20[rank*2+2])
        else:
            text(obs.index[-1]+ datetime.timedelta(hours=0.1), yposlab[rank], column, fontsize=14, color=tableau20[rank+2])
        
    for y in range(10, upbound, 10):  
        axhline(y, lw=0.5, color="black", alpha=0.3)
    
    text(obs.index[len(obs)/2], upbound+2, ''.join(["Soil Moisture and Precipitation at ",catch]) , fontsize=17, ha="center") 
    text(obs.index[0]- datetime.timedelta(hours=0.1),32,'< Soil Moisture\n',rotation='vertical',horizontalalignment='right',verticalalignment='bottom',fontsize=12,alpha=0.7)
    text(obs.index[0]- datetime.timedelta(hours=0.1),32,'> Precipitation',rotation='vertical',horizontalalignment='right',verticalalignment='bottom',fontsize=12,color=tableau20[0],alpha=0.7)
    text(obs.index[-1]+ datetime.timedelta(hours=0.1),10,descriptor,rotation='vertical',verticalalignment='bottom',fontsize=12,alpha=0.7)
    text(obs.index[-1]+ datetime.timedelta(hours=0.1),7.55,'Hr',verticalalignment='bottom',fontsize=12,    alpha=0.7)
    text(obs.index[-1]+ datetime.timedelta(hours=0.1),5.6,'Date',verticalalignment='bottom',fontsize=14,    alpha=0.7)
        
    ax2 = ax.twinx()
    if cumprec:
        obs['Preccum']=np.cumsum(obs['Prec'].values)
        obs=obs.rename(columns={'Prec': 'Prec_m','Preccum': 'Prec'})

    precstp=np.around(obs['Prec'].max()/4.,decimals=-int(np.floor(np.log10(obs['Prec'].max()/4.))))
    ax2.set_ylim((-obs['Prec'].max()*3.,1))
    ax2.set_yticks(-np.arange(0,obs['Prec'].max(), precstp)[::-1])
    ax2.set_yticklabels([str(x) + "mm" for x in np.arange(0.,obs['Prec'].max(), precstp)/precscale][::-1],    fontsize=14,color=tableau20[0])
    ax2.plot(obs.index, -obs['Prec'], color=tableau20[0])
    ax2.fill_between(obs.index, -obs['Prec'], 0., color=tableau20[1])
    ax2.yaxis.grid(True)
    
    if saving!=False:
        savefig(saving, bbox_inches="tight")


    

In [5]:

    

mlab=['0.03 m','0.1 m','0.2 m','0.3 m','0.4 m']
mlabcols=['BF_3','BF_10','BF_20','BF_30','BF_40']
#hydroplot(obs,mlab,mlabcols,(10,6),'weiher_moist.pdf',catch='Weiherbach')
hydroplot(obs,mlab,mlabcols,(10,6),catch='Weiherbach',align=True)


    

In [6]:

    

reftime=pd.to_datetime('21:35 27/06/94')
idx=np.argmin([np.abs(obs.index[x]-reftime) for x in range(len(obs))])
#obs.iloc[idx:idx+20]


    

In [7]:

    

import scipy as sp
import matplotlib.pyplot as plt
import os, sys


    

In [8]:

    

#connect SimPEG Flow
pathdir='../'
lib_path = os.path.abspath(pathdir)
sys.path.append(lib_path)
pathdir='../../../../code/simpegflow/'
lib_path = os.path.abspath(pathdir)
sys.path.append(lib_path)

import vG_conv as vG
from SimPEG import *
from simpegFLOW import Richards


    

    




Warning: upgrade your scipy to 0.13.0

In [9]:

    

#connect to echoRD
import run_echoRD as rE
#connect and load project
[dr,mc,mcp,pdyn,cinf,vG]=rE.loadconnect(pathdir='../',mcinif='mcini_specht4y',oldvers=True)
mc.advectref='obs'
[mc,particles,npart,precTS]=rE.pickup_echoRD(mc, mcp, dr, 'specht4y.pickle')


    

    




updated mc from specht4y.pickle

In [10]:

    

#define grid resolution
grid_res=-mc.mgrid.vertfac.values #[m]
#setup mesh
cells=mc.mgrid.vertgrid.values
M = Mesh.TensorMesh([np.ones(cells)*grid_res])
M.setCellGradBC([['neumann','dirichlet']])


    

    
Out[10]:





[['neumann', 'dirichlet']]

In [11]:

    

#assign soil parameters
params = {}
#access soil parameters from simpegFlow for initialisation
vgParams = Richards.Empirical.VanGenuchtenParams()
for param in vgParams.clayLoam:
    params[param] =  vgParams.clayLoam[param]*np.ones(cells)

params['Ks'] = np.log(params['Ks']) #correct the given Ks values


    

In [12]:

    

#update params with our soil (homogeneous for now)
params['Ks'] = np.log(mc.soilmatrix.ks)[4].repeat(len(params['Ks']))
params['alpha'] = mc.soilmatrix.alpha[4].repeat(len(params['alpha']))
params['n'] = mc.soilmatrix.n[4].repeat(len(params['n']))
params['theta_r'] = mc.soilmatrix.tr[4].repeat(len(params['theta_r']))
params['theta_s'] = mc.soilmatrix.ts[4].repeat(len(params['theta_s']))


    

In [13]:

    

#pre-event moisture
#reftime=pd.to_datetime('23:55 27/06/94')
idx=np.argmin([np.abs(obs.index[x]-reftime) for x in range(len(obs))])
obs.iloc[idx]


    

    
Out[13]:





BF_3     24.070000
BF_10    18.100000
BF_20    18.780001
BF_30    22.770000
BF_40    27.730000
Prec      0.000000
Name: 1994-06-27 21:35:00, dtype: float64

In [14]:

    

#impose measurement to setup
theta=mc.zgrid[:,1]*np.nan
probloc=[-0.03,-0.1,-0.2,-0.3,-0.4]
for i in range(len(probloc)):
    j=np.argmin(np.abs(mc.zgrid[:,1]-probloc[i]))
    theta[j]=obs.iloc[idx,i]
theta[0]=obs.iloc[idx,0]
theta[-1]=obs.iloc[idx,i]
theta=pd.DataFrame(theta).interpolate().values[:,0]
theta=theta/100.


    

In [15]:

    

#psi profile
E = Richards.Empirical.VanGenuchten(M, **params)
h=vG.psi_theta(theta, E.theta_s, E.theta_r,E.alpha,E.n)
mtrue = (params['Ks'])


    

In [16]:

    

#build boundary conditions - take over values
bcBTx = np.array([h[-1],h[0]])
def bc(time,top_psi):
    bcBT = bcBTx
    return bcBT


    

In [17]:

    

ths_part=mc.part_sizefac
npart=np.floor(mc.part_sizefac*vG.thst_theta(theta,mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])).values.astype(int)
#dgrid=np.append(-np.diff(mc.zgrid[:,1]),-np.diff(mc.zgrid[0:2,1]))
#npart=np.floor(dgrid*ths_part*vG.thst_theta(theta,mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])).values.astype(int)
particles=pd.DataFrame(np.zeros(int(np.sum(npart))*9).reshape(int(np.sum(npart)),9),columns=['lat', 'z', 'conc', 'temp', 'age', 'flag', 'fastlane', 'advect','lastZ'])
particles['cell']=pd.Series(np.zeros(int(np.sum(npart)),dtype=int), index=particles.index)
k=0

cells=len(npart)

for i in np.arange(len(npart)):
    particles.cell[k:(k+npart[i])]=int(i)
    particles.lat[k:(k+npart[i])]=(np.random.rand(npart[i]))*mc.mgrid.latfac.values
    particles.z[k:(k+npart[i])]=(i+np.random.rand(npart[i]))*mc.mgrid.vertfac.values
    k+=npart[i]


    

In [18]:

    

def gridupdate_thS1D(cellid,mc,pdyn):
    '''Calculates thetaS from particle density
    '''
    cells=np.append(range(mc.mgrid.cells),cellid)
    cells[cells<0]=0
    cells[cells>len(mc.zgrid[:,1])]=len(mc.zgrid[:,1])
    thetaS=np.floor(100.*np.bincount(particles.cell).astype(float)/mc.part_sizefac)
    thetaS[thetaS>100.]=100.
    return [thetaS.astype(np.int),npart]

def part_diffusion_1D(particles,npart,thS,mc,vG,pdyn,dt,uffink_corr=True):
    '''Calculate Diffusive Particle Movement
       Based on state in grid use diffusivity as foundation of 2D random walk.
       Project step and check for boundary conditions and further restrictions.
       Update particle positions.
    '''
    N=len(particles.z) #number of particles

    # 1D Random Walk function with additional correction term for
    # non-static diffusion after Uffink 1990 p.15 & p.24ff and Kitanidis 1994
    xi=np.random.rand(N)*2.-1.

    u=mc.ku[thS,mc.soilgrid[:,4]-1]/mc.theta[thS,mc.soilgrid[:,4]-1]
    D=mc.D[thS,mc.soilgrid[:,4]-1]*mc.soilmatrix.ts[mc.soilgrid[:,4]-1].values
 
    vert_sproj=(dt*u[particles.cell.values.astype(np.int)] + (xi*((2*D[particles.cell.values.astype(np.int)]*dt)**0.5)))

    if (uffink_corr==True):
        #Itô Scheme after Uffink 1990 and Kitanidis 1994 for vertical step
        #modified Stratonovich Scheme after Kitanidis 1994 for lateral step
        dx=vert_sproj
        
        # project step and updated state
        # new positions
        z_proj=particles.z.values-vert_sproj
        [lat_proj,z_proj,nodrain]=pdyn.boundcheck(particles.lat,z_proj,mc)
        [thSx,npartx]=gridupdate_thS1D(pdyn.cellgrid(particles.lat,z_proj,mc).astype(np.int64),mc,pdyn) 
        cell_proj=pdyn.cellgrid(particles.lat,z_proj,mc).astype(np.int64)
        
        u_proj=mc.ku[thSx,mc.soilgrid[:,4]-1]/mc.theta[thSx,mc.soilgrid[:,4]-1]
        D_proj=mc.D[thSx,mc.soilgrid[:,4]-1]*mc.soilmatrix.ts[mc.soilgrid[:,4]-1].values
    
        corrD=np.abs(D_proj[particles.cell.values.astype(np.int)]-D[particles.cell.values.astype(np.int)])/dx
        corrD[dx==0.]=0.
        D_mean=np.sqrt(D_proj[particles.cell.values.astype(np.int)]*D[particles.cell.values.astype(np.int)])
        corru=np.sqrt(u[particles.cell.values.astype(np.int)]*u_proj[particles.cell.values.astype(np.int)])

        vert_sproj=((corru-corrD)*dt + (xi*((2*D[particles.cell.values.astype(np.int)]*dt)**0.5)))


    # new positions
    z_new=particles.z.values-vert_sproj

    [particles.lat,particles.z,nodrain]=pdyn.boundcheck(particles.lat,z_proj,mc)
    particles['cell']=pdyn.cellgrid(particles.lat,particles.z,mc).astype(np.int64)

    # saturation check
    [thS,npart]=gridupdate_thS1D(pdyn.cellgrid(particles.lat,z_new,mc).astype(np.int64),mc,pdyn) #DEBUG: externalise smooth parameter

    if any(-nodrain):
        particles.loc[-nodrain,'flag']=len(mc.maccols)+1

    return [particles,thS,npart]

def part_advect_1D(particles,dt):
    idx=particles['flag']!=0
    if any(idx):
        particles.z.loc[idx]=particles.z.loc[idx]+particles.advect.loc[idx]*dt
        particles.cell.loc[idx]=pdyn.cellgrid(particles.lat.loc[idx],particles.z.loc[idx],mc).astype(np.int64).values
    return particles

def infilt_1Dobs(ti,precip,mc,pdyn,dt):
    T=np.array(9)
    # get timestep in prec time series
    prec_id=np.argmin([np.abs(precip.index[x]-ti) for x in range(len(precip))])
    if precip.iloc[prec_id]>0:
        prec_avail=int(np.round(precip.iloc[prec_id]/600.*dt*mc.part_sizefac/(-mc.mgrid.vertfac.values*1000.)))
        #avail. water particles
        prec_c=np.array([23.])
    else:
        prec_avail=0
        prec_c=np.array([0.,0.])
    
    if prec_avail>0:
        particles_infilt=pd.DataFrame(np.zeros(prec_avail*10).reshape(prec_avail,10),columns=['lat', 'z', 'conc', 'temp', 'age', 'flag', 'fastlane', 'advect','lastZ','cell'])
        # place particles at surface and redistribute later according to ponding
        particles_infilt.z=-0.00001
        particles_infilt.lat=np.random.rand(prec_avail)*mc.mgrid.latfac.values
        particles_infilt.conc=prec_c[0]
        particles_infilt.temp=T
        particles_infilt.age=ti
        particles_infilt.flag=1
        particles_infilt.fastlane=np.random.randint(len(mc.t_cdf_fast.T), size=prec_avail)
        particles_infilt.cell=0
        particles_infilt.advect=pdyn.assignadvect(prec_avail,mc,particles_infilt.fastlane.values,False)
    else:
        particles_infilt=pd.DataFrame([])

    return particles_infilt


    

In [19]:

    

def CAOSpy_run1D_adv(particles,npart,thS,leftover,drained,tstart,tstop,precTS,mc,pdyn,cinf):
    timenow=tstart
    precparts=0
    #loop through time
    while timenow < tstop:
        #define dt as Courant/Neumann criterion
        dt_D=(mc.mgrid.vertfac.values[0])**2 / (np.amax(mc.D[np.amax(thS),:]))
        dt_ku=-mc.mgrid.vertfac.values[0]/np.amax(mc.ku[np.amax(thS),:])
        dt=np.amin([dt_D,dt_ku])
        
        #INFILT
        p_inf=infilt_1Dobs(timenow,precTS,mc,pdyn,dt)
        particles=pd.concat([particles,p_inf])
        #psi=vG.psi_thst(thS/100.,mc.soilmatrix.alpha[mc.soilgrid[:,1]-1],mc.soilmatrix.n[mc.soilgrid[:,1]-1])
        #DIFFUSION
        [particles,thS,npart]=part_diffusion_1D(particles,npart,thS,mc,vG,pdyn,dt,uffink_corr=True)
        
        #ADVECTION
        particles=part_advect_1D(particles,dt)
        
        #drained particles
        drained=drained.append(particles[particles.flag==len(mc.maccols)+1])
        particles=particles[particles.flag!=len(mc.maccols)+1]
        pondparts=(particles.z<0.)
        leftover=np.count_nonzero(-pondparts)
        [particles.lat,particles.z,nodrain]=pdyn.boundcheck(particles.lat,particles.z,mc)
        particles['cell']=pdyn.cellgrid(particles.lat,particles.z,mc).astype(np.int64)
        [thS,npart]=gridupdate_thS1D(particles.cell,mc,pdyn)
        timenow=timenow+dt
        precparts+=len(p_inf)

    return(particles,npart,thS,leftover,drained,timenow,precparts)


    

In [20]:

    

#some plotting functions

def oneDplot(particles,obsx,theta_r,theta_re,dt,sigma,runname,ti,i,mc,saving=False,t=1,store=False,fsize=(8, 5),xlimset=[0.15,0.3,2],ad_diff=False,noplot=False):
    #plot 1D profile
    import matplotlib.gridspec as gridspec
    from scipy.ndimage.filters import gaussian_filter1d
    [thS,npart]=gridupdate_thS1D(particles.cell,mc,pdyn)
    theta_p=vG.theta_thst(thS/100., mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])
    thpx=gaussian_filter1d(theta_p,sigma)
    if ad_diff:
        [thSdiff,npartdiff]=gridupdate_thS1D(particles.cell[particles.flag==0],mc,pdyn)
        theta_pdiff=vG.theta_thst(thSdiff/100., mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])
        thpxdiff=gaussian_filter1d(theta_pdiff,sigma)
    n=len(particles)
    obs_id=np.argmin([np.abs(obsx.index[x]-ti) for x in range(len(obsx))])
    probloc=[-0.03,-0.1,-0.2,-0.3,-0.4]
    if -noplot:
        fig=plt.figure(figsize=fsize)
        gs = gridspec.GridSpec(1, 2, width_ratios=[2,1])
        subplot(gs[0])
        plot(thpx,mc.zgrid[:,1],label='Particle')
        if ad_diff:
            plot(thpxdiff,mc.zgrid[:,1],label='Particle_diffusive')
        plot(theta_r,mc.zgrid[:,1],label='Rich SimpegFlow')
        plot(theta_re,mc.zgrid[:,1],label='Rich Euler')
        plot(obsx.iloc[obs_id]/100.,probloc,'.',label='Observation')
        legend(loc=4)
        #text(0.35, -0.4, ''.join(['t=',str(int(t)),'m']), fontsize=12)
        #text(0.3, -0.5, ''.join(['particles: ',str(n)]), fontsize=12)
        xlim(xlimset[:2])
        xticks(np.arange(xlimset[0],xlimset[1],xlimset[2]))
        xlabel('theta [m3/m3]')
        ylabel('depth [m]')
        #title(''.join(['Model and Obs @ ',str(int(ti)),'s']))
        title(''.join(['Model and Observation\nTime: ',str(int(t)),'min']))
        #title(''.join(['echoRD1D @ ',str(int(ti)),'s']))
        ax1=subplot(gs[1])
        zi=np.arange(-0.0,mc.soildepth-0.01,-0.01)
        oldp=np.bincount(np.round(np.append(-particles.z[particles.flag==0].values,-zi)*100.).astype(int))-1
        allp=np.bincount(np.round(np.append(-particles.z[particles.flag<2].values,-zi)*100.).astype(int))-1
        plot(gaussian_filter1d(oldp,sigma),zi,label='old')
        plot(gaussian_filter1d(allp[0:len(oldp)],sigma),zi,label='all')
        plot(gaussian_filter1d(allp[0:len(oldp)]-oldp,sigma),zi,label='new')
        a=np.ceil(n/1000.)*12.
        xlim([0,a])
        xticks(np.linspace(0,a,4))
        xlabel('particles')
        legend(loc=4)
        #title(''.join(['Max Peclet=',str(np.round(Pe,2))]))
        title(''.join(['total:\n',str(n)]))
        ax1.yaxis.tick_right()
    
        if saving:
            plt.savefig(''.join(['./results/',runname,str(i).zfill(3),'.pdf']))
        
        plt.close(fig)
        
    if store:
        idz=[0,10,20,30,40]
        if ad_diff:
            return [obsx.values[obs_id]/100., thpx[idz],theta_re.values[idz],theta_r[idz],thpxdiff[idz]]
        else:
            return [obsx.values[obs_id]/100., thpx[idz],theta_re.values[idz],theta_r[idz]]
    
        
def oneDplot2(particles,obsx,theta_r,theta_re,dt,sigma,runname,ti,i,mc,saving=False,t=1,store=False,fsize=(8, 5),xlimset=[0.15,0.3,2],ad_diff=False):
    #plot 1D profile
    import matplotlib.gridspec as gridspec
    from scipy.ndimage.filters import gaussian_filter1d
    [thS,npart]=gridupdate_thS1D(particles.cell,mc,pdyn)
    theta_p=vG.theta_thst(thS/100., mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])
    thpx=gaussian_filter1d(theta_p,sigma)
    if ad_diff:
        [thSdiff,npartdiff]=gridupdate_thS1D(particles.cell[particles.flag==0],mc,pdyn)
        theta_pdiff=vG.theta_thst(thSdiff/100., mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])
        thpxdiff=gaussian_filter1d(theta_pdiff,sigma)
    n=len(particles)
    obs_id=np.argmin([np.abs(obsx.index[x]-ti) for x in range(len(obsx))])
    probloc=[-0.03,-0.1,-0.2,-0.3,-0.4]
    
    fig=plt.figure(figsize=fsize)
    gs = gridspec.GridSpec(1, 4, width_ratios=[2,1,0.6,0.6])
    subplot(gs[0])
    plot(thpx,mc.zgrid[:,1],label='Particle')
    if ad_diff:
        plot(thpxdiff,mc.zgrid[:,1],label='Particle_diffusive')
    plot(theta_r,mc.zgrid[:,1],label='Rich SimpegFlow')
    plot(theta_re,mc.zgrid[:,1],label='Rich Euler')
    plot(obsx.iloc[obs_id]/100.,probloc,'.',label='Observation')
    legend(loc=4)
    #text(0.35, -0.4, ''.join(['t=',str(int(t)),'m']), fontsize=12)
    #text(0.3, -0.5, ''.join(['particles: ',str(n)]), fontsize=12)
    xlim(xlimset[:2])
    xticks(np.arange(xlimset[0],xlimset[1],xlimset[2]))
    xlabel('theta [m3/m3]')
    ylabel('depth [m]')
    #title(''.join(['Model and Obs @ ',str(int(ti)),'s']))
    title(''.join(['Model and Observation\nTime: ',str(int(t)),'min']))
    #title(''.join(['echoRD1D @ ',str(int(ti)),'s']))
    ax1=subplot(gs[1])
    zi=np.arange(-0.0,mc.soildepth-0.01,-0.01)
    oldp=np.bincount(np.round(np.append(-particles.z[particles.flag==0].values,-zi)*100.).astype(int))-1
    allp=np.bincount(np.round(np.append(-particles.z[particles.flag<2].values,-zi)*100.).astype(int))-1
    plot(gaussian_filter1d(oldp,sigma),zi,label='old')
    plot(gaussian_filter1d(allp[0:len(oldp)],sigma),zi,label='all')
    plot(gaussian_filter1d(allp[0:len(oldp)]-oldp,sigma),zi,label='new')
    a=np.ceil(n/1000.)*12.
    xlim([0,a])
    xticks(np.linspace(0,a,4))
    xlabel('particles')
    legend(loc=4)
    #title(''.join(['Max Peclet=',str(np.round(Pe,2))]))
    title(''.join(['total:\n',str(n)]))
    ax1.get_yaxis().set_visible(False)    
    
    ax2=subplot(gs[2])
    #a=(np.abs(particles.groupby(['cell'], sort=True).max().lastZ.values-particles.groupby(['cell'], sort=True).min().z.values)/dt)/np.sqrt(mc.D[np.amax(thS),4])
    #e=particles.groupby(['cell'], sort=True).mean().lastD.values*mc.mgrid.vertfac.values/mc.D[np.amax(thS),4]
    a=(np.abs(mc.mgrid.vertfac.values)*np.abs(particles.groupby(['cell'], sort=True).mean().lastZ.values-particles.groupby(['cell'], sort=True).mean().z.values)/dt)/mc.D[np.amax(thS),4]
    e=(np.abs(mc.mgrid.vertfac.values)*np.abs(particles.groupby(['cell'], sort=True).max().lastZ.values-particles.groupby(['cell'], sort=True).min().z.values)/dt)/mc.D[np.amax(thS),4]
    plot(a,mc.zgrid[:,1],label='diff')
    #plot(e,mc.zgrid[:,1],alpha=0.5, color='b')
    #fill_betweenx(mc.zgrid[:,1], a, e,alpha=0.2)
    #plot(u,mc.zgrid[:,1],color='adv')
    #errorbar(a,mc.zgrid[:,1],xerr=a-e, ecolor='lightblue')
    xlim([-0.05,np.ceil(np.amax([np.amax(a),0.05])*100.)/100.])
    xticks(np.linspace(0.,np.ceil(np.amax([np.amax(a),0.05])*100.)/100.,2))
    ax2.get_yaxis().set_visible(False)
    #xscale('log')
    xlabel('u(i,z)/D(z)')
    title('Peclet\n(mean)')
    
    ax3=subplot(gs[3])
    plot(e,mc.zgrid[:,1])
    xlim([-0.05,np.ceil(np.amax([np.amax(e),0.05])*100.)/100.])
    xticks(np.linspace(0.,np.ceil(np.amax([np.amax(e),0.05])*100.)/100.,2))
    #xscale('log')
    xlabel('u(i,z)/D(z)')
    title('Peclet\n(max)')
    ax3.yaxis.tick_right()
    
    if saving:
        plt.savefig(''.join(['./results/',runname,str(i).zfill(3),'.pdf']))
        plt.close(fig)
    if store:
        idz=[0,10,20,30,40]
        if ad_diff:
            return [obsx.values[obs_id]/100., thpx[idz],theta_re.values[idz],theta_r[idz],thpxdiff[idz]]
        else:
            return [obsx.values[obs_id]/100., thpx[idz],theta_re.values[idz],theta_r[idz]]
        
def oneDplotX(particles,obsx,theta_r,theta_re,dt,sigma,runname,ti,i,mc,saving=False,t=1,store=False,fsize=(8, 5),xlimset=[0.15,0.3,2],leg=False,ad_diff=False,relative=False):
    #plot 1D profile
    import matplotlib.gridspec as gridspec
    from scipy.ndimage.filters import gaussian_filter1d
    
    # These are the "Tableau 20" colors as RGB.  
    tableau20 = [(31, 119, 180), (174, 199, 232), (255, 127, 14), (255, 187, 120),  
             (44, 160, 44), (152, 223, 138), (214, 39, 40), (255, 152, 150),  
             (148, 103, 189), (197, 176, 213), (140, 86, 75), (196, 156, 148),  
             (227, 119, 194), (247, 182, 210), (127, 127, 127), (199, 199, 199),  
             (188, 189, 34), (219, 219, 141), (23, 190, 207), (158, 218, 229)]  
    # Scale the RGB values to the [0, 1] range, which is the format matplotlib accepts.  
    for k in range(len(tableau20)):  
        r, g, b = tableau20[k]  
        tableau20[k] = (r / 255., g / 255., b / 255.)  
    
    [thS,npart]=gridupdate_thS1D(particles.cell,mc,pdyn)
    theta_p=vG.theta_thst(thS/100., mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])
    thpx=gaussian_filter1d(theta_p,sigma)
    if ad_diff:
        [thSdiff,npartdiff]=gridupdate_thS1D(particles.cell[particles.flag==0],mc,pdyn)
        theta_pdiff=vG.theta_thst(thSdiff/100., mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])
        thpxdiff=gaussian_filter1d(theta_pdiff,sigma)
    n=len(particles)
    obs_id=np.argmin([np.abs(obsx.index[x]-ti) for x in range(len(obsx))])
    probloc=[-0.03,-0.1,-0.2,-0.3,-0.4]
    
    fig=plt.figure(figsize=fsize)
    gs = gridspec.GridSpec(1, 2, width_ratios=[4,1])
    
    #marginal Y
    ax1=fig.add_subplot(gs[1])
    zi=np.arange(-0.0,mc.soildepth-0.01,-0.01)
    oldp=np.bincount(np.round(np.append(-particles.z[particles.flag==0].values,-zi)*100.).astype(int))-1
    allp=np.bincount(np.round(np.append(-particles.z[particles.flag<2].values,-zi)*100.).astype(int))-1
    if relative:
        oldp/=np.sum(oldp)
        allp/=np.sum(allp)
    advect_dummy=gaussian_filter1d(10.*(allp[0:len(oldp)]-oldp),sigma)
    old_dummy=gaussian_filter1d(oldp,sigma)
    all_dummy=gaussian_filter1d(allp[0:len(oldp)],sigma)
    ax1.fill_betweenx(zi,0.,all_dummy, color='b',alpha=0.15)
    ax1.fill_betweenx(zi,0.,old_dummy, color='g',alpha=0.15)
    ax1.fill_betweenx(zi,0.,advect_dummy, color='r',alpha=0.3)
    ax1.plot(advect_dummy,zi,'r-',label='new particles')
    ax1.plot(all_dummy,zi,'b-',label='all particles')
    ax1.plot(old_dummy,zi,'g-',label='old particles')
    
    handles1, labels1 = ax1.get_legend_handles_labels() 
    ax1.set_ylim((-1.,0.))
    ax1.set_yticks([])
    ax1.set_xticks([])
    ax1.spines["bottom"].set_visible(False)  
    ax1.spines["right"].set_visible(False)  
    ax1.spines["top"].set_visible(False)  
    ax1.spines["left"].set_visible(False)
    
    #main plot
    ax1=fig.add_subplot(gs[0])
    ax1.plot(thpx,mc.zgrid[:,1],color=tableau20[0],lw=2,label='Particle')
    if ad_diff:
        ax1.plot(thpxdiff,mc.zgrid[:,1],'--',color=tableau20[1],lw=2,label='Particle_diffusive')
    ax1.plot(theta_r,mc.zgrid[:,1],color=tableau20[2],lw=2,label='Rich SimpegFlow')
    ax1.plot(theta_re,mc.zgrid[:,1],color=tableau20[4],lw=2,label='Rich Euler')
    ax1.plot(obsx.iloc[obs_id]/100.,probloc,'o',color=tableau20[8],label='Observation')
    
    ax1.set_xlim(xlimset[:2])
    ax1.set_xticks(np.arange(xlimset[0],xlimset[1],xlimset[2]))
    ax1.set_xticklabels(np.arange(xlimset[0],xlimset[1],xlimset[2]),fontsize=14)
    ax1.set_yticklabels([-1.0,'',-0.6,-0.4,-0.2,0.0],fontsize=14)
    ax1.text(0.12,-0.75,'depth [m]',rotation='vertical',fontsize=14)
    if leg:
        ax1.legend(loc=4,frameon=False)
    ax1.set_xlabel('theta [m3/m3]',fontsize=14)
    ax1.text(0.17,-0.95,''.join([str(int(i)),'min']),fontsize=19)
        
    if saving:
        plt.savefig(''.join(['./results/',runname,str(i).zfill(3),'.pdf']))
        plt.close(fig)


    

In [21]:

    

# a simple progress bar
import sys, time
try:
    from IPython.display import clear_output
    have_ipython = True
except ImportError:
    have_ipython = False

class ProgressBar:
    def __init__(self, iterations):
        self.iterations = iterations
        self.prog_bar = '[]'
        self.fill_char = '*'
        self.width = 40
        self.__update_amount(0)
        if have_ipython:
            self.animate = self.animate_ipython
        else:
            self.animate = self.animate_noipython

    def animate_ipython(self, iter):
        print '\r', self,
        sys.stdout.flush()
        self.update_iteration(iter + 1)

    def update_iteration(self, elapsed_iter):
        self.__update_amount((elapsed_iter / float(self.iterations)) * 100.0)
        self.prog_bar += '  %d of %s complete' % (elapsed_iter, self.iterations)

    def __update_amount(self, new_amount):
        percent_done = int(round((new_amount / 100.0) * 100.0))
        all_full = self.width - 2
        num_hashes = int(round((percent_done / 100.0) * all_full))
        self.prog_bar = '[' + self.fill_char * num_hashes + ' ' * (all_full - num_hashes) + ']'
        pct_place = (len(self.prog_bar) // 2) - len(str(percent_done))
        pct_string = '%d%%' % percent_done
        self.prog_bar = self.prog_bar[0:pct_place] + \
            (pct_string + self.prog_bar[pct_place + len(pct_string):])

    def __str__(self):
        return str(self.prog_bar)


    

In [22]:

    

#euler richards solver
def darcy(psi,z,k):
    phi=psi+z
    dz=np.diff(z)
    Q=np.sqrt(k[:-1]*k[1:])*np.diff(phi)/dz
    return Q

def euler(psi_l,c_l,q,qb_u,qb_l,dt,z):
    dz=np.diff(z)
    #psinew=psi_l[:-1]+dt/(c_l[:-1]*dz*dz)*np.append(qb_u,np.diff(q))
    psinew=psi_l[:-1]+dt/(c_l[:-1]*dz*dz)*np.diff(np.append((qb_u+q[0])/2.,q))
    psilow=psi_l[-1]+dt/(c_l[-1]*dz[-1])*qb_l
    return np.append(psinew,psilow)

def richards(t_end,psi,mc):
    time=0.
    soil=mc.soilgrid[:,1]-1
    dzmn=mc.mgrid.latfac.values
    while time<t_end:
        k=vG.ku_psi(psi, mc.soilmatrix.ks[soil], mc.soilmatrix.alpha[soil], mc.soilmatrix.n[soil]).values
        c=vG.c_psi(psi, mc.soilmatrix.ts[soil], mc.soilmatrix.tr[soil], mc.soilmatrix.alpha[soil], mc.soilmatrix.n[soil]).values
        q=darcy(psi,mc.zgrid[:,1],k)
        dt=np.amin([0.1,0.05*dzmn/np.amax(np.abs(q))])
        #predictor
        psinew=euler(psi,c,q,0.,0.,dt*0.5,mc.zgrid[:,1])

        k=vG.ku_psi(psinew, mc.soilmatrix.ks[soil], mc.soilmatrix.alpha[soil], mc.soilmatrix.n[soil]).values
        c=vG.c_psi(psinew, mc.soilmatrix.ts[soil], mc.soilmatrix.tr[soil], mc.soilmatrix.alpha[soil], mc.soilmatrix.n[soil]).values
        q=darcy(psinew,mc.zgrid[:,1],k)
    
        #corrector
        psinew=euler(psinew,c,q,0.,0.,dt,mc.zgrid[:,1])
        psi=psinew
        time=time+dt
    
    return psi


    

In [23]:

    

#prepare for 1D
mc.mgrid.latgrid=1
mc.mgrid.width=mc.mgrid.latfac*40.
drained=pd.DataFrame(np.array([]))
leftover=0
runname='echoRD_1D_xobs_semidry10X'


    

In [24]:

    

[thS,npart]=gridupdate_thS1D(particles.cell,mc,pdyn)
endtime=pd.to_datetime('06:00 28/06/94')
t_end=(endtime-reftime).seconds
output=60.
dummy=np.floor(t_end/output)
dummi=dummy.astype(int)
leftover=0
drained=pd.DataFrame(np.array([]))
timenow=0.
p = ProgressBar(dummi)
mc.prects=True
ad_diff=False


    

In [25]:

    

#prepare precip time series for echoRD
idy=np.argmin([np.abs(obs.index[x]-endtime) for x in range(len(obs))])
precTSx=obs.Prec.iloc[idx:idy]/100.
precTSx.index=[np.abs(precTSx.index[x]-reftime).seconds for x in range(len(precTSx))]


    

In [26]:

    

#prepare obs time series
obsx=obs.iloc[idx:idy,0:5]
obsx.index=[np.abs(obsx.index[x]-reftime).seconds for x in range(len(obsx))]


    

In [27]:

    

soil=mc.soilgrid[:,1]-1


    

In [28]:

    

#output=1000.
#test simpegFlow
h = vG.psi_theta(theta,mc.soilmatrix.ts[soil],mc.soilmatrix.tr[soil],mc.soilmatrix.alpha[soil],mc.soilmatrix.n[soil]).values
def bc(time,top_psi):
    return np.array([h[0],h[-1]])
prob = Richards.RichardsProblem(M, mapping=E, timeSteps=[(output, 1)], boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed')
Hs = prob.fields(mtrue)
theta_r=E.theta(Hs[1],mtrue)

#test eulerRichards
psi_re = richards(output,h,mc)
theta_re = vG.theta_psi(psi_re, mc.soilmatrix.ts[soil], mc.soilmatrix.tr[soil], mc.soilmatrix.alpha[soil], mc.soilmatrix.n[soil])
oneDplot(particles,obsx,theta_r,theta_re,0.,1,runname,timenow,i,mc,False,output*i/60.)